# This script generates figures A1 and A2 and tables 1, 2 and A1-A17.

rm(list = ls())


library(here)

# Define relative path to save in the script's directory
file_path <- here()


load(here('crossnational_data.rda'))

library(zoo)
library(dplyr)
library(magrittr)
library(kableExtra)
library(countrycode)
library(panelView)
library(ggplot2)
library(fixest)
library(modelsummary)



# Interpolating some of the missing values between years here.
crossnational_data <- crossnational_data %>% group_by(Country) %>% mutate(e_wbgi_cce_imp = na.approx(e_wbgi_cce, na.rm=FALSE)) %>% ungroup()  


# Create GDP percent change
crossnational_data <- crossnational_data %>% dplyr::group_by(Country) %>% mutate(gdp_change_pp = ((e_gdp*100) / dplyr::lag(e_gdp))-100)

# Create a descriptive table (Table A1) ----

important_covariates <- crossnational_data %>% 
  ungroup() %>% 
  dplyr::select(Satis, SupDem, v2x_jucon, v2xlg_legcon, v2xcl_rol, gdp_change_pp, gini_disp, e_peaveduc, e_wbgi_cce_imp, v2x_clphy)

variable_labels <- c(
  Satis = "Satisfaction with democracy",
  SupDem = 'Support for democracy',
  v2x_jucon = "Judicial constraints",
  v2xlg_legcon = 'Legislative constraints',
  v2xcl_rol = 'Individual liberty',
  gdp_change_pp = 'GDP change (percent)',
  gini_disp = 'Gini index (disposable)',
  e_peaveduc = 'Years of education (ave.)',
  e_wbgi_cce_imp = 'Control of corruption',
  v2x_clphy = 'Physical violence index'
)

# Calculate the descriptive statistics
descriptive_table <- important_covariates %>%
  summarise(across(everything(), list(
    N = ~sum(!is.na(.)),
    Mean = ~round(mean(., na.rm = TRUE), 2),
    Sd = ~round(sd(., na.rm = TRUE), 2),
    Min = ~round(min(., na.rm = TRUE), 2),
    Max = ~round(max(., na.rm = TRUE), 2)
  ), .names = "{.col}_{.fn}")) %>%  
  # Ensure the column names are correctly formatted
  tidyr::pivot_longer(
    cols = everything(), 
    names_to = c("variable", "stat"), 
    names_pattern = "(.*)_(.*)"  # Regex to split at the last underscore
  ) %>% 
  tidyr::pivot_wider(names_from = "stat", values_from = "value") %>% 
  mutate(variable = recode(variable, !!!variable_labels)) %>% 
  rename(Variable = variable)


# Generate plain text table
descriptive_table <- kable( descriptive_table, format = "pipe")

# Write the descriptive table:
writeLines(descriptive_table, here('table_A1.txt'))


#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Data coverage plot ----

crossnational_data$country_code <- countrycode(crossnational_data$Country, origin = "country.name",  destination = "iso3c")

crossnational_data$country_code <- ifelse(crossnational_data$Country == 'Kosovo', 'XKX', crossnational_data$country_code)

crossnational_data$coverage <- ifelse(is.na(crossnational_data$Satis)==0, 1, 0)

# Sort countries alphabetically
crossnational_data <- crossnational_data[order(crossnational_data$Country),]

# Creating two different figures to show data coverage: 

# Find the midpoint of unique countries
countries <- unique(crossnational_data$Country)
midpoint <- ceiling(length(countries) / 2)

# Create subsets
crossnational_data_1 <- crossnational_data[crossnational_data$Country %in% countries[1:midpoint], ]

crossnational_data_2 <- crossnational_data[crossnational_data$Country %in% countries[(midpoint + 1):length(countries)], ]


first_plot <- panelview(Satis ~ coverage, data = crossnational_data_1, index = c("country_code", "Year")) + theme_bw() + theme(legend.position = 'bottom', axis.text.x = element_text(size = 8, angle = 90, vjust = 0.5, hjust = 1), axis.text.y = element_text(size = 6), legend.title = element_blank()) + labs(title = '', y = 'Country Code')

ggsave(first_plot, file = here('figure_A1.pdf'), height = 7, width = 6)

second_plot <- panelview(Satis ~ coverage, data = crossnational_data_2, index = c("country_code", "Year")) + theme_bw() + theme(legend.position = 'bottom', axis.text.x = element_text(size = 8, angle = 90, vjust = 0.5, hjust = 1), axis.text.y = element_text(size = 6), legend.title = element_blank()) + labs(title = '', y = 'Country Code')


ggsave(second_plot, file = here('figure_A2.pdf'), height = 7, width = 6)


#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Dependent variable: Democracy Satisfaction ----

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

## Independent variable: Judicial constraints ----

m1 <- feols(Satis~ v2x_jucon | Country, data=crossnational_data, vcov = ~Country)

m2 <- feols(Satis~ v2x_jucon | Country + Year, data=crossnational_data, vcov = ~Country)


m3 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp | Country, data=crossnational_data, vcov = ~Country)

m4 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp | Country+Year, data=crossnational_data, vcov = ~Country)


m5 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country, data=crossnational_data, vcov = ~Country)

m6 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=crossnational_data, vcov = ~Country)


m7 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=crossnational_data, vcov = ~Country)

m8 <- feols(Satis~ v2x_jucon+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=crossnational_data, vcov = ~Country)


m9 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=crossnational_data, vcov = ~Country)

m10 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country + Year, data=crossnational_data, vcov = ~Country)

## substantive interpretation:
abs(sd(crossnational_data$v2x_jucon, na.rm = TRUE)*coef(m10)['v2x_jucon']) / sd(crossnational_data$Satis, na.rm = TRUE)

abs(sd(crossnational_data$e_wbgi_cce_imp, na.rm = TRUE)*coef(m10)['e_wbgi_cce_imp']) / sd(crossnational_data$Satis, na.rm = TRUE)

# models for appendix:
jucon_appendix1 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy+e_miinflat | Country, data=crossnational_data, vcov = ~Country)

jucon_appendix2 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy+e_miinflat | Country + Year, data=crossnational_data, vcov = ~Country)

jucon_appendix3 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy+e_miinflat+e_total_oil_income_pc | Country, data=crossnational_data, vcov = ~Country)

jucon_appendix4 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy+e_miinflat+e_total_oil_income_pc | Country + Year, data=crossnational_data, vcov = ~Country)



modellist <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)



model_rows <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`, ~`9`, ~`10`,
                      'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes',
                      'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')

attr(model_rows, 'position') <- c(13, 14)


# Generate plain text output
model_table <- modelsummary( modellist,  stars_note = TRUE,  stars = c('*' = .1, '**' = .05, '***' = 0.01),  gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', 
coef_map = c(
    'v2x_jucon' = 'Judicial constraints', 
    'gdp_change_pp' = 'GDP change (percent)', 
    'gini_disp' = 'Gini index (disposable)', 
    'e_peaveduc' = 'Years of education (ave.)', 
    'e_wbgi_cce_imp' = 'Control of corruption', 
    'e_miinflat' = 'Inflation', 
    'v2x_clphy' = 'Physical violence index', 
    'e_total_oil_income_pc' = 'Petroleum production per capita'), 
  add_rows = model_rows, 
  notes = 'Standard errors are clustered at the country level',  
  output = here("table_1.txt"))


## Independent variable: Legislative constraints ----

m1 <- feols(Satis~ v2xlg_legcon | Country, data=crossnational_data, vcov = ~Country)

m2 <- feols(Satis~ v2xlg_legcon | Country + Year, data=crossnational_data, vcov = ~Country)


m3 <- feols(Satis~ v2xlg_legcon+gdp_change_pp+gini_disp | Country, data=crossnational_data, vcov = ~Country)

m4 <- feols(Satis~ v2xlg_legcon+gdp_change_pp+gini_disp | Country+Year, data=crossnational_data, vcov = ~Country)


m5 <- feols(Satis~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc | Country, data=crossnational_data, vcov = ~Country)

m6 <- feols(Satis~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=crossnational_data, vcov = ~Country)


m7 <- feols(Satis~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=crossnational_data, vcov = ~Country)

m8 <- feols(Satis~ v2xlg_legcon+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=crossnational_data, vcov = ~Country)


m9 <- feols(Satis~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=crossnational_data, vcov = ~Country)

m10 <- feols(Satis~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country+Year, data=crossnational_data, vcov = ~Country)


## Models for Appendix:
legcon_appendix1 <- feols(Satis~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy+e_miinflat | Country, data=crossnational_data, vcov = ~Country)

legcon_appendix2 <- feols(Satis~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy+e_miinflat | Country + Year, data=crossnational_data, vcov = ~Country)

legcon_appendix3 <- feols(Satis~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy+e_miinflat+e_total_oil_income_pc | Country, data=crossnational_data, vcov = ~Country)

legcon_appendix4 <- feols(Satis~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy+e_miinflat+e_total_oil_income_pc | Country + Year, data=crossnational_data, vcov = ~Country)



modellist <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)


modelsummary(modellist, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2xlg_legcon' = 'Legislative constraints', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows, notes = 'Standard errors are clustered at the country level',  output = here("table_A3.txt"))

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@


## Independent variable: judicial and legislative contraints ----

modellist_appendix <- list('1' = jucon_appendix1, '2' = jucon_appendix2, '3' = jucon_appendix3, '4' = jucon_appendix4, '5' = legcon_appendix1, '6' = legcon_appendix2, '7' = legcon_appendix3, '8' = legcon_appendix4)

model_rows_appendix <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`,
                               'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes',
                               'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')


modelsummary(modellist_appendix, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2x_jucon' = 'Judicial constraints', 'v2xlg_legcon' = 'Legislative contraints', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows_appendix, notes = 'Standard errors are clustered at the country level',  output = here("table_A2.txt"))

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

## Independent variable: Ind liberties ----

m1 <- feols(Satis~ v2xcl_rol | Country, data=crossnational_data, vcov = ~Country)

m2 <- feols(Satis~ v2xcl_rol | Country + Year, data=crossnational_data, vcov = ~Country)

m3 <- feols(Satis~ v2xcl_rol+gdp_change_pp+gini_disp | Country, data=crossnational_data, vcov = ~Country)

m4 <- feols(Satis~ v2xcl_rol+gdp_change_pp+gini_disp | Country+Year, data=crossnational_data, vcov = ~Country)


m5 <- feols(Satis~ v2xcl_rol+gdp_change_pp+gini_disp + e_peaveduc | Country, data=crossnational_data, vcov = ~Country)

m6 <- feols(Satis~ v2xcl_rol+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=crossnational_data, vcov = ~Country)


m7 <- feols(Satis~ v2xcl_rol+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=crossnational_data, vcov = ~Country)

m8 <- feols(Satis~ v2xcl_rol+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=crossnational_data, vcov = ~Country)


m9 <- feols(Satis~ v2xcl_rol+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=crossnational_data, vcov = ~Country)

m10 <- feols(Satis~ v2xcl_rol+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country+Year, data=crossnational_data, vcov = ~Country)


modellist <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)


modelsummary(modellist, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2xcl_rol' = 'Individual liberty', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows, notes = 'Standard errors are clustered at the country level',  output = here("table_2.txt")) 

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Dependent variable: Democracy Mood ----

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@


## Independent variable: Judicial constraints ----

m1 <- feols(SupDem~ v2x_jucon | Country, data=crossnational_data, vcov = ~Country)

m2 <- feols(SupDem~ v2x_jucon | Country + Year, data=crossnational_data, vcov = ~Country)


m3 <- feols(SupDem~ v2x_jucon+gdp_change_pp+gini_disp | Country, data=crossnational_data, vcov = ~Country)

m4 <- feols(SupDem~ v2x_jucon+gdp_change_pp+gini_disp | Country+Year, data=crossnational_data, vcov = ~Country)


m5 <- feols(SupDem~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country, data=crossnational_data, vcov = ~Country)

m6 <- feols(SupDem~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=crossnational_data, vcov = ~Country)


m7 <- feols(SupDem~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=crossnational_data, vcov = ~Country)

m8 <- feols(SupDem~ v2x_jucon+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=crossnational_data, vcov = ~Country)


m9 <- feols(SupDem~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=crossnational_data, vcov = ~Country)

m10 <- feols(SupDem~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country + Year, data=crossnational_data, vcov = ~Country)


modellist <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)



model_rows <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`, ~`9`, ~`10`,
                      'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes',
                      'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')

attr(model_rows, 'position') <- c(13, 14)


# Generate plain text output
model_table <- modelsummary( modellist,  stars_note = TRUE,  stars = c('*' = .1, '**' = .05, '***' = 0.01),  gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', 
                             coef_map = c(
                               'v2x_jucon' = 'Judicial constraints', 
                               'gdp_change_pp' = 'GDP change (percent)', 
                               'gini_disp' = 'Gini index (disposable)', 
                               'e_peaveduc' = 'Years of education (ave.)', 
                               'e_wbgi_cce_imp' = 'Control of corruption', 
                               'e_miinflat' = 'Inflation', 
                               'v2x_clphy' = 'Physical violence index', 
                               'e_total_oil_income_pc' = 'Petroleum production per capita'), 
                             add_rows = model_rows, 
                             notes = 'Standard errors are clustered at the country level',  
                             output = here("table_A6.txt"))


## Independent variable: Legislative constraints ----

m1 <- feols(SupDem~ v2xlg_legcon | Country, data=crossnational_data, vcov = ~Country)

m2 <- feols(SupDem~ v2xlg_legcon | Country + Year, data=crossnational_data, vcov = ~Country)


m3 <- feols(SupDem~ v2xlg_legcon+gdp_change_pp+gini_disp | Country, data=crossnational_data, vcov = ~Country)

m4 <- feols(SupDem~ v2xlg_legcon+gdp_change_pp+gini_disp | Country+Year, data=crossnational_data, vcov = ~Country)


m5 <- feols(SupDem~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc | Country, data=crossnational_data, vcov = ~Country)

m6 <- feols(SupDem~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=crossnational_data, vcov = ~Country)


m7 <- feols(SupDem~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=crossnational_data, vcov = ~Country)

m8 <- feols(SupDem~ v2xlg_legcon+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=crossnational_data, vcov = ~Country)


m9 <- feols(SupDem~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=crossnational_data, vcov = ~Country)

m10 <- feols(SupDem~ v2xlg_legcon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country+Year, data=crossnational_data, vcov = ~Country)



modellist <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)


modelsummary(modellist, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2xlg_legcon' = 'Legislative constraints', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows, notes = 'Standard errors are clustered at the country level',  output = here("table_A7.txt"))


#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

## Independent variable: Ind liberties ----

m1 <- feols(SupDem~ v2xcl_rol | Country, data=crossnational_data, vcov = ~Country)

m2 <- feols(SupDem~ v2xcl_rol | Country + Year, data=crossnational_data, vcov = ~Country)

m3 <- feols(SupDem~ v2xcl_rol+gdp_change_pp+gini_disp | Country, data=crossnational_data, vcov = ~Country)

m4 <- feols(SupDem~ v2xcl_rol+gdp_change_pp+gini_disp | Country+Year, data=crossnational_data, vcov = ~Country)


m5 <- feols(SupDem~ v2xcl_rol+gdp_change_pp+gini_disp + e_peaveduc | Country, data=crossnational_data, vcov = ~Country)

m6 <- feols(SupDem~ v2xcl_rol+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=crossnational_data, vcov = ~Country)


m7 <- feols(SupDem~ v2xcl_rol+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=crossnational_data, vcov = ~Country)

m8 <- feols(SupDem~ v2xcl_rol+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=crossnational_data, vcov = ~Country)


m9 <- feols(SupDem~ v2xcl_rol+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=crossnational_data, vcov = ~Country)

m10 <- feols(SupDem~ v2xcl_rol+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country+Year, data=crossnational_data, vcov = ~Country)


modellist <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)


modelsummary(modellist, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2xcl_rol' = 'Individual liberty', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows, notes = 'Standard errors are clustered at the country level',  output = here("table_A8.txt")) 

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Alternative independent variables: Democracy Satisfaction ----

## Horizontal Accountability ----

m1 <- feols(Satis~ v2x_horacc | Country, data=crossnational_data, vcov = ~Country)

m2 <- feols(Satis~ v2x_horacc | Country + Year, data=crossnational_data, vcov = ~Country)


m3 <- feols(Satis~ v2x_horacc+gdp_change_pp+gini_disp | Country, data=crossnational_data, vcov = ~Country)

m4 <- feols(Satis~ v2x_horacc+gdp_change_pp+gini_disp | Country+Year, data=crossnational_data, vcov = ~Country)


m5 <- feols(Satis~ v2x_horacc+gdp_change_pp+gini_disp + e_peaveduc | Country, data=crossnational_data, vcov = ~Country)

m6 <- feols(Satis~ v2x_horacc+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=crossnational_data, vcov = ~Country)


m7 <- feols(Satis~ v2x_horacc+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=crossnational_data, vcov = ~Country)

m8 <- feols(Satis~ v2x_horacc+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=crossnational_data, vcov = ~Country)


m9 <- feols(Satis~ v2x_horacc+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=crossnational_data, vcov = ~Country)

m10 <- feols(Satis~ v2x_horacc+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country + Year, data=crossnational_data, vcov = ~Country)


modellist_horacc <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)

model_rows_horacc <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`, ~`9`, ~`10`,
                             'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes',
                             'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')


attr(model_rows_horacc, 'position') <- c(13, 14)


modelsummary(modellist_horacc, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2x_horacc' = 'Horizontal accountability index', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows_horacc, notes = 'Standard errors are clustered at the country level', output = here("table_A4.txt")) 

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

## Vertical Accountability ----

m1 <- feols(Satis~ v2x_veracc | Country, data=crossnational_data, vcov = ~Country)

m2 <- feols(Satis~ v2x_veracc | Country + Year, data=crossnational_data, vcov = ~Country)

m3 <- feols(Satis~ v2x_veracc+gdp_change_pp+gini_disp | Country, data=crossnational_data, vcov = ~Country)

m4 <- feols(Satis~ v2x_veracc+gdp_change_pp+gini_disp | Country+Year, data=crossnational_data, vcov = ~Country)


m5 <- feols(Satis~ v2x_veracc+gdp_change_pp+gini_disp + e_peaveduc | Country, data=crossnational_data, vcov = ~Country)

m6 <- feols(Satis~ v2x_veracc+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=crossnational_data, vcov = ~Country)


m7 <- feols(Satis~ v2x_veracc+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=crossnational_data, vcov = ~Country)

m8 <- feols(Satis~ v2x_veracc+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=crossnational_data, vcov = ~Country)


m9 <- feols(Satis~ v2x_veracc+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=crossnational_data, vcov = ~Country)

m10 <- feols(Satis~ v2x_veracc+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country + Year, data=crossnational_data, vcov = ~Country)


modellist_veracc <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)


model_rows_veracc <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`, ~`9`, ~`10`,
                             'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes',
                             'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')


attr(model_rows_veracc, 'position') <- c(13, 14)

modelsummary(modellist_veracc, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2x_veracc' = 'Vertical accountability index', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows_veracc, notes = 'Standard errors are clustered at the country level',  output = here("table_A5.txt")) 

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Alternative independent variables:  Mood ----

## Horizontal Accountability ----

m1 <- feols(SupDem~ v2x_horacc | Country, data=crossnational_data, vcov = ~Country)

m2 <- feols(SupDem~ v2x_horacc | Country + Year, data=crossnational_data, vcov = ~Country)

m3 <- feols(SupDem~ v2x_horacc+gdp_change_pp+gini_disp | Country, data=crossnational_data, vcov = ~Country)

m4 <- feols(SupDem~ v2x_horacc+gdp_change_pp+gini_disp | Country+Year, data=crossnational_data, vcov = ~Country)


m5 <- feols(SupDem~ v2x_horacc+gdp_change_pp+gini_disp + e_peaveduc | Country, data=crossnational_data, vcov = ~Country)

m6 <- feols(SupDem~ v2x_horacc+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=crossnational_data, vcov = ~Country)


m7 <- feols(SupDem~ v2x_horacc+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=crossnational_data, vcov = ~Country)

m8 <- feols(SupDem~ v2x_horacc+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=crossnational_data, vcov = ~Country)


m9 <- feols(SupDem~ v2x_horacc+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=crossnational_data, vcov = ~Country)

m10 <- feols(SupDem~ v2x_horacc+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country + Year, data=crossnational_data, vcov = ~Country)


modellist_horacc <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)


model_rows_horacc <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`, ~`9`, ~`10`,
                             'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes',
                             'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')


attr(model_rows_horacc, 'position') <- c(13, 14)


modelsummary(modellist_horacc, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2x_horacc' = 'Horizontal accountability index', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows_horacc, notes = 'Standard errors are clustered at the country level',  output = here("table_A9.txt")) 

## Vertical accountability ----

m1 <- feols(SupDem~ v2x_veracc | Country, data=crossnational_data, vcov = ~Country)

m2 <- feols(SupDem~ v2x_veracc | Country + Year, data=crossnational_data, vcov = ~Country)


m3 <- feols(SupDem~ v2x_veracc+gdp_change_pp+gini_disp | Country, data=crossnational_data, vcov = ~Country)

m4 <- feols(SupDem~ v2x_veracc+gdp_change_pp+gini_disp | Country+Year, data=crossnational_data, vcov = ~Country)


m5 <- feols(SupDem~ v2x_veracc+gdp_change_pp+gini_disp + e_peaveduc | Country, data=crossnational_data, vcov = ~Country)

m6 <- feols(SupDem~ v2x_veracc+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=crossnational_data, vcov = ~Country)


m7 <- feols(SupDem~ v2x_veracc+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=crossnational_data, vcov = ~Country)

m8 <- feols(SupDem~ v2x_veracc+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=crossnational_data, vcov = ~Country)


m9 <- feols(SupDem~ v2x_veracc+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=crossnational_data, vcov = ~Country)

m10 <- feols(SupDem~ v2x_veracc+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country + Year, data=crossnational_data, vcov = ~Country)


modellist_veracc <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)


model_rows_veracc <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`, ~`9`, ~`10`,
                             'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes',
                             'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')


attr(model_rows_veracc, 'position') <- c(13, 14)


modelsummary(modellist_veracc, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2x_veracc' = 'Vertical accountability index', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows_veracc, notes = 'Standard errors are clustered at the country level',  output = here("table_A10.txt")) 

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Main results heterogeneity ----

crossnational_data$continent <- ifelse(crossnational_data$e_regiongeo %in% c(1, 2, 3, 4), 'Europe', ifelse(crossnational_data$e_regiongeo %in% c(5, 6, 7, 8, 9), 'Africa', ifelse(crossnational_data$e_regiongeo %in% c(10, 11, 12, 13, 14), 'Asia', ifelse(crossnational_data$e_regiongeo == 15, 'Ocenia', ifelse(crossnational_data$e_regiongeo %in% c(16, 17, 18, 19), 'America', NA)))))

europe_sample <- crossnational_data %>% filter(continent == 'Europe')

america_sample <- crossnational_data %>% filter(continent == 'America')

asia_sample <- crossnational_data %>% filter(continent == 'Asia')

africa_sample <- crossnational_data %>% filter(continent == 'Africa')


#### America sample ----

m1 <- feols(Satis~ v2x_jucon | Country, data=america_sample, vcov = ~Country)

m2 <- feols(Satis~ v2x_jucon | Country + Year, data=america_sample, vcov = ~Country)

m3 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp | Country, data=america_sample, vcov = ~Country)

m4 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp | Country+Year, data=america_sample, vcov = ~Country)


m5 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country, data=america_sample, vcov = ~Country)

m6 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=america_sample, vcov = ~Country)


m7 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=america_sample, vcov = ~Country)

m8 <- feols(Satis~ v2x_jucon+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=america_sample, vcov = ~Country)


m9 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=america_sample, vcov = ~Country)

m10 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country + Year, data=america_sample, vcov = ~Country)


modellist <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)



model_rows <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`, ~`9`, ~`10`,
                      'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes',
                      'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')

attr(model_rows, 'position') <- c(13, 14)



modelsummary(modellist, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2x_jucon' = 'Judicial constraints', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows, notes = 'Standard errors are clustered at the country level', output = here("table_A11.txt")) 

#### Europe sample ----


m1 <- feols(Satis~ v2x_jucon | Country, data=europe_sample, vcov = ~Country)

m2 <- feols(Satis~ v2x_jucon | Country + Year, data=europe_sample, vcov = ~Country)


m3 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp | Country, data=europe_sample, vcov = ~Country)

m4 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp | Country+Year, data=europe_sample, vcov = ~Country)


m5 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country, data=europe_sample, vcov = ~Country)

m6 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=europe_sample, vcov = ~Country)


m7 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=europe_sample, vcov = ~Country)

m8 <- feols(Satis~ v2x_jucon+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=europe_sample, vcov = ~Country)


m9 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=europe_sample, vcov = ~Country)

m10 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country + Year, data=europe_sample, vcov = ~Country)

modellist <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)



model_rows <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`, ~`9`, ~`10`,
                      'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes',
                      'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')

attr(model_rows, 'position') <- c(13, 14)


modelsummary(modellist, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2x_jucon' = 'Judicial constraints', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows, notes = 'Standard errors are clustered at the country level',  output = here("table_A14.txt")) 

#### Asia sample ----


m1 <- feols(Satis~ v2x_jucon | Country, data=asia_sample, vcov = ~Country)

m2 <- feols(Satis~ v2x_jucon | Country + Year, data=asia_sample, vcov = ~Country)


m3 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp | Country, data=asia_sample, vcov = ~Country)

m4 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp | Country+Year, data=asia_sample, vcov = ~Country)


m5 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country, data=asia_sample, vcov = ~Country)

m6 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=asia_sample, vcov = ~Country)


m7 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=asia_sample, vcov = ~Country)

m8 <- feols(Satis~ v2x_jucon+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=asia_sample, vcov = ~Country)


m9 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=asia_sample, vcov = ~Country)

m10 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country + Year, data=asia_sample, vcov = ~Country)



modellist <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)


model_rows <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`, ~`9`, ~`10`,
                      'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes',
                      'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')

attr(model_rows, 'position') <- c(13, 14)

modelsummary(modellist, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2x_jucon' = 'Judicial constraints', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows, notes = 'Standard errors are clustered at the country level',  output = here("table_A12.txt")) 

#### Africa sample ----


m1 <- feols(Satis~ v2x_jucon | Country, data=africa_sample, vcov = ~Country)

m2 <- feols(Satis~ v2x_jucon | Country + Year, data=africa_sample, vcov = ~Country)


m3 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp | Country, data=africa_sample, vcov = ~Country)

m4 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp | Country+Year, data=africa_sample, vcov = ~Country)


m5 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country, data=africa_sample, vcov = ~Country)

m6 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=africa_sample, vcov = ~Country)


m7 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=africa_sample, vcov = ~Country)

m8 <- feols(Satis~ v2x_jucon+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=africa_sample, vcov = ~Country)


m9 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=africa_sample, vcov = ~Country)

m10 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country + Year, data=africa_sample, vcov = ~Country)


modellist <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)


model_rows <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`, ~`9`, ~`10`,
                      'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes',
                      'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')

attr(model_rows, 'position') <- c(13, 14)


modelsummary(modellist, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2x_jucon' = 'Judicial constraints', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows, notes = 'Standard errors are clustered at the country level',  output = here("table_A13.txt")) 

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

## Liberal vs. electoral democracies ----

crossnational_data <- crossnational_data %>% arrange(Country, Year) %>% group_by(Country) %>% 
  mutate(v2x_jucon_lag = lag(v2x_jucon, n = 1),
         gdp_change_pp_lag = lag(gdp_change_pp, n = 1), 
         gini_disp_lag = lag(gini_disp, n = 1),
         e_peaveduc_lag = lag(e_peaveduc, n = 1),
         e_wbgi_cce_imp_lag = lag(e_wbgi_cce_imp, n = 1),
         v2x_clphy_lag = lag(v2x_clphy, n = 1)) %>% ungroup()  

lib_democ <- crossnational_data %>% filter(v2x_regime %in% c(3))

elec_democ <- crossnational_data %>% filter(v2x_regime %in% c(2))

democracies <- crossnational_data %>% filter(v2x_regime %in% c(2, 3))

lib_1 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=lib_democ, vcov = ~Country)

lib_2 <- feols(Satis~ v2x_jucon+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=lib_democ, vcov = ~Country)

lib_3 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=lib_democ, vcov = ~Country)

lib_4 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country + Year, data=lib_democ, vcov = ~Country)

elec_1 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=elec_democ, vcov = ~Country)

elec_2 <- feols(Satis~ v2x_jucon+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=elec_democ, vcov = ~Country)


elec_3 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=elec_democ, vcov = ~Country)

elec_4 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country + Year, data=elec_democ, vcov = ~Country)

modellist <- list('1' = lib_1, '2' = lib_2, '3' = lib_3, '4' = lib_4, '5' = elec_1, '6' = elec_2, '7' = elec_3, '8' = elec_4)


model_rows <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`,
                      'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 
                      'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')

attr(model_rows, 'position') <- c(13, 14)

modelsummary(modellist, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2x_jucon' = 'Judicial constraints', 'gdp_change_pp' = 'GDP change (percent)', 'gini_disp' = 'Gini index (disposable)', 'e_peaveduc' = 'Years of education (ave.)', 'e_wbgi_cce_imp' = 'Control of corruption', 'e_miinflat' = 'Inflation', 'v2x_clphy' = 'Physical violence index', 'e_total_oil_income_pc' = 'Petroleum production per capita'), add_rows = model_rows, notes = 'Standard errors are clustered at the country level',   output = here("table_A16.txt"))


## Lagged indp vars:

lib_1 <- feols(Satis ~ v2x_jucon_lag +gdp_change_pp_lag +gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag | Country, data=lib_democ, vcov = ~Country)

lib_2 <- feols(Satis ~ v2x_jucon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag +e_wbgi_cce_imp_lag | Country + Year, data=lib_democ, vcov = ~Country)

lib_3 <- feols(Satis ~ v2x_jucon_lag +gdp_change_pp_lag +gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag + v2x_clphy_lag | Country, data=lib_democ, vcov = ~Country)

lib_4 <- feols(Satis~ v2x_jucon_lag +gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag + v2x_clphy_lag | Country + Year, data=lib_democ, vcov = ~Country)


elec_1 <- feols(Satis~ v2x_jucon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag | Country, data=elec_democ, vcov = ~Country)

elec_2 <- feols(Satis~ v2x_jucon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag | Country + Year, data=elec_democ, vcov = ~Country)


elec_3 <- feols(Satis~ v2x_jucon_lag + gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag + v2x_clphy_lag | Country, data=elec_democ, vcov = ~Country)

elec_4 <- feols(Satis~ v2x_jucon_lag +gdp_change_pp_lag + gini_disp_lag + e_peaveduc_lag + e_wbgi_cce_imp_lag + v2x_clphy_lag | Country + Year, data=elec_democ, vcov = ~Country)

modellist <- list('1' = lib_1, '2' = lib_2, '3' = lib_3, '4' = lib_4, '5' = elec_1, '6' = elec_2, '7' = elec_3, '8' = elec_4)


model_rows <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`,
                      'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 
                      'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')

attr(model_rows, 'position') <- c(13, 14)

modelsummary(modellist, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', coef_map = c('v2x_jucon_lag' = 'Judicial constraints$_{t-1}$', 'gdp_change_pp_lag' = 'GDP change (percent)$_{t-1}$', 'gini_disp_lag' = 'Gini index (disposable)$_{t-1}$', 'e_peaveduc_lag' = 'Years of education (ave.)$_{t-1}$', 'e_wbgi_cce_imp_lag' = 'Control of corruption$_{t-1}$', 'e_miinflat_lag' = 'Inflation$_{t-1}$', 'v2x_clphy_lag' = 'Physical violence index$_{t-1}$', 'e_total_oil_income_pc_lag' = 'Petroleum production per capita$_{t-1}$'), add_rows = model_rows,  output = here("table_A17.txt"), escape = FALSE)


#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

## Democracy sample only ----

m1 <- feols(Satis~ v2x_jucon | Country, data=democracies, vcov = ~Country)

m2 <- feols(Satis~ v2x_jucon | Country + Year, data=democracies, vcov = ~Country)

m3 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp | Country, data=democracies, vcov = ~Country)

m4 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp | Country+Year, data=democracies, vcov = ~Country)


m5 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country, data=democracies, vcov = ~Country)

m6 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc | Country + Year, data=democracies, vcov = ~Country)


m7 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country, data=democracies, vcov = ~Country)

m8 <- feols(Satis~ v2x_jucon+ gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp | Country + Year, data=democracies, vcov = ~Country)


m9 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country, data=democracies, vcov = ~Country)

m10 <- feols(Satis~ v2x_jucon+gdp_change_pp+gini_disp + e_peaveduc+e_wbgi_cce_imp+v2x_clphy | Country + Year, data=democracies, vcov = ~Country)


modellist <- list('1' = m1, '2' = m2, '3' = m3, '4' = m4, '5' = m5, '6' = m6, '7' = m7, '8' = m8, '9' = m9, '10' = m10)



model_rows <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`, ~`9`, ~`10`,
                      'Country FE', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes', 'Yes',
                      'Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes')

attr(model_rows, 'position') <- c(13, 14)


# Generate plain text output
model_table <- modelsummary( modellist,  stars_note = TRUE,  stars = c('*' = .1, '**' = .05, '***' = 0.01),  gof_omit = 'DF|Deviance|AIC|BIC|F|Log.Lik.|R2 Pseudo|R2 Within|Std.Errors|RMSE', 
                             coef_map = c(
                               'v2x_jucon' = 'Judicial constraints', 
                               'gdp_change_pp' = 'GDP change (percent)', 
                               'gini_disp' = 'Gini index (disposable)', 
                               'e_peaveduc' = 'Years of education (ave.)', 
                               'e_wbgi_cce_imp' = 'Control of corruption', 
                               'e_miinflat' = 'Inflation', 
                               'v2x_clphy' = 'Physical violence index', 
                               'e_total_oil_income_pc' = 'Petroleum production per capita'), 
                             add_rows = model_rows, 
                             notes = 'Standard errors are clustered at the country level',  
                             output = here("table_A15.txt"))

